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Abstract 

The models that are based of fractional derivatives should be highlighted among 
promising new models to describe turbulent fluid flows. In the present work, a steady- 
state flow in a duct is considered under the condition that the turbulent diffusion is 
governed by a fractional power of the Laplace operator. To study numerically flows in 
rectangular channels, finite-difference approximations are employed. For approximate 
solving the corresponding boundary value problem, the iterative method of conjugate 
gradients is used. At each iteration, the problem with a fractional power of the grid 
Laplace operator is solved. Predictions of turbulent flows in ducts at different Reynolds 
numbers are presented via mean velocity fields. 

Keywords: Turbulent flow, fluid flow in ducts, fractional power of the Laplace 
operator, finite-difference problem, iterative method of conjugate gradients 


1. Introduction 

To model continuum mechanics phenomena, different models of turbulence are em¬ 
ployed (see, e.g., Jl]|2] among others). In terms of practical use, emphasis is on simple 
mathematical models of turbulence, which, on the one hand, are not much more com¬ 
plex in comparison with models for laminar flows, and on the other hand, reproduce 
the basic features of turbulent regimes of liquid and gas flows. 

Nowadays, non-local applied mathematical models based on using fractional deriva¬ 
tives in time and space are actively discussed |[3][4j|5]. Many models in applied physics, 
biology, hydrology and finance, involve both sub-diffusion (fractional in time) and 
supper-diffusion (fractional in space) operators. Supper-diffusion problems are treated 
as evolutionary problems with a fractional power of an elliptic operator. 

Such an anomalous diffusion model is used in |j6) to describe turbulent flows. In the 
work Q, a turbulent diffusion in the Reynolds equations for the mean velocity is gov- 
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erned by the fractional Laplacian. The development of this approach is hindered by the 
lack of simple and robust numerical algorithms for solving boundary value problems 
for equations with fractional powers. In the best case (see, e.g., ||8]), investigations are 
restricted to simple one-dimensional in space models. 

For solving problems with fractional powers of elliptic operators, we can apply 
finite volume and finite element methods oriented to using arbitrary domains and ir¬ 
regular computational grids mm- The numerical implementation involves the matrix 
function-vector multiplication. For such problems, different approaches CD are avail¬ 
able. Problems of using Krylov subspace methods with the Lanczos approximation 
when solving systems of linear equations associated with the fractional elliptic equa¬ 
tions are discussed in D2. a comparative analysis of the contour integral method, the 
extended Krylov subspace method, and the preassigned poles and interpolation nodes 
method for solving space-fractional reaction-diffusion equations is presented in ED- 
The simplest variant is associated with the explicit construction of the solution using 
the known eigenvalues and eigenfunctions of the elliptic operator with diagonalization 
of the corresponding matrix msa. Unfortunately, all these approaches demonstrates 
too high computational complexity for multidimensional problems. 

We have proposed QU a computational algorithm for solving an equation with 
fractional powers of elliptic operators on the basis of a transition to a pseudo-parabolic 
equation. For the auxiliary Cauchy problem, the standard two-level schemes are ap¬ 
plied. The computational algorithm is simple for practical use, robust, and applicable 
to solving a wide class of problems. A small number of time steps is required to find 
a solution. This computational algorithm for solving equations with fractional powers 
of operators is promising when considering transient problems. 

To implement numerically a space-fractional model of turbulent fluid flow, we must 
take into account a multi-term structure of the problem operator. Namely, here one 
term is the standard elliptic operator (normal diffusion), whereas the second term is 
a fractional power of an elliptic operator (anomalous diffusion). For solving such 
non-classical problems, it seems natural to apply iterative methods with an appropriate 
choice of preconditioners mm. 

In this paper, for predicting a steady-state turbulent flow in a duct, we apply a model 
with a turbulent space-fractional diffusion. To solve numerically this problem with the 
multi-term diffusion, we employ the iterative method of conjugate gradients, where the 
problem with normal diffusion is solved at each iteration to construct a preconditioner. 
For solving the problem with the fractional Laplacian, a pseudo-parabolic equation is 
used. The paper is organized as follows. In Section 2, a mathematical model with the 
fractional Laplacian is introduced to describe a turbulent flow in a rectangular duct. 
The discrete problem and computational algorithm are discussed in Section 3. Section 
4 presents an analysis of the impact of the basic parameters of the problem on numerical 
results obtained using the developed model. 

2. A space-fractional model of turbulent fluid flow 

Motion of an incompressible fluid is governed by the Navier-Stokes equations: 
dv 1 

— + v ■ Vv + -'Vp - vAv = 0, (1) 

at p 


2 


V • V = 0. 


(2) 


Here p is the density, p denotes the pressure, v stands for the velocity vector, and v is 
the fluid viscosity. 

To obtain the Reynolds equations for turbulent flows [T|[2], velocity and pres¬ 
sure v, p are decomposed into the sum of the mean flow components v, p and fluc¬ 
tuating components V, ~p. Substituting this decomposition into (JTJ, 0, we arrive at 
the Reynolds equations written in the following coordinate-wise representation (v = 
(vi,v 2 ,v 3 )): 


dvj _ dvi l dp _ d 

-h V;~ -1 -VAV; + —— VjVj = 0, 

dt J dxj p dxi dxi 


(3) 



(4) 


A RANS model of turbulence is defined by a particular formulation for the Reynolds 
stress tensor pvjv}. 

For the space-fractional model, we have 


—Wj = ?(-A )°Vi. (5) 

Here the coefficient is treated as the eddy (turbulent) diffusivity. In the work E), 
some arguments are given in favor of setting the power a equal to 1/3. In view of in¬ 
equations 0-0 may be written in the form similar to 0, ([ 2 ]), i.e., 

dv 1 

— +v- Vv+ -Vp- VAV+^(-A) ff V = 0, ( 6 ) 

dt p 


V ■ v = 0. (7) 

Let us consider a steady-state stabilized in the longitudinal direction flow in rect¬ 
angular channels (x = ( X\,X 2 )): 


O = {jc | x = (x\ , x 2 ) 0 < Xj < dj, i = 1,2). 

Let x 3 be the longitudinal coordinate and assume that v = (0,0, u). Then from 0, 0, 
we obtain the following equation for the longitudinal component of the velocity: 

- vAm + %(-A) a u - x, x e fl, (8) 


where 


, s 1 dp 

P=p(x 3 ), X = - -!—■ 

p ax 3 


The equation 0 is supplemented with homogeneous Dirichlet boundary conditions: 

u(x) = 0, x e <9Q, (9) 


which corresponds to the no-slip condition on rigid walls. 
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For the normalization of equation ([8|, as the reference values, we employ the chan¬ 
nel height d 2 and the velocity scale 


M 0 


4 


-x- 


For the dimensionless velocity n, using for the dimensionless quantities the same nota¬ 
tion as for the dimensional ones, we obtain 


Am + /j(-A) a u =1, x e Q, 


( 10 ) 


where 

Q = {x | x - (x\,x 2 ) | 0 < xi < d, 0 < X 2 < 1), 

f ,2(1-0-) 

F = - d 2 ■ 
v z 

Thus, the boundary value problem (j9j), ( |T()| has three governing parameters, namely, 
a , ju and d. 

3. Computational algorithm 

To solve the steady-state problem (|9j, we introduce a uniform grid in the 
domain O: 


to - [x \ x - (xi,X 2 ), Xk = ikhk, h - 0,1, ...,Nk, N\h\ = d, A^ 2^2 = 1), 

with Td — a) U dot, where lo is the set of interior points and da> is the set of boundary 
grid points. For grid functions y(x) such that vfx) = 0, x t tu, we define the Hilbert 
space // = L 2 (tu), where the scalar product and the norm are given as follows: 

(y,w) = 'Y j y(x)w{x)hih 2 , ||y|| = (y,y ) 1/2 . 

xeco 

For the discrete Laplace operator A, we introduce the additive representation 

2 

A = ^ A k , x e u>, (11) 

k=l 

where Ak, k = 1,2 are associated with the corresponding differential operator of the 
second derivative in one direction. 

For all grid points except adjacent to the boundary, the grid operator A\ can be 
written as 


My = - (>’(■* i 


- /Ji,/i 2 ) - 2y(x) +y(xi - h u h 2 )). 


x 6 cj, x, + 0.5/zi, xi + d - 0.5hi. 
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In the points that are adjacent to the boundary, the approximation is constructed taking 
into account the boundary condition |[9|: 


A\y = -4r(y(jci + h\Ji 2 ) - 2y(jc)), jtew, x\ = 0.5/ii, 

h\ 

A\y =—(2y(x) - y(x\ - h\,h 2 )), x e a>, xi=d-0.5h\. 

h- 

Similarly we construct the grid operator ,4 2 . For the above grid operators, we have (see, 

e.g„ HU ED) 

A k =A\> 6 k E , 4 = 4 sin 2 * = 1,2. 

h k 2N k 


where E is the identity operator. Because of this, the discrete Laplace operator (11 1 is 
self-adjoint and positive definite in H: 


2 

A - A* > 6E, d = 

k=l 


( 12 ) 


It approximates the differential Laplace operator with the truncation error <9(|/z| 2 ), 
\h\ 2 = h\ + h\. 

To handle the fractional power of the grid operator A , let us consider the eigenvalue 
problem 

A^Ptn — 

which has the well-known analytical solution. We have 


6 = di < A 2 < ... < A m , M = (N\ - 1 )(N 2 - 1), 


where eigenfunctions ip m , ||i/3,„|| = 1, m = 1,2,..., M form a basis in H. Therefore 

M 

y = ^jA,‘Pm)Vm- (13) 

m= 1 

For the fractional power of the operator A, we have 

M 

Ay— 'y ^ (y, ip m )A m ip m . 

m= 1 

Using the above approximations, we arrive from (|9j, ( [I()| at the discrete problem 

Ay + fiA a y = 1. (14) 


In our particular case with using uniform meshes in a rectangle, the solution of equation 
( [T4| can be constructed explicitly via the known eigenvalues and eigenfunctions. For 
the solution represented in the form of (p~3]>, we obtain 


0 , Vm) 


(1 , ^Pm ) 
Am + l-lA m 


M. 
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We are interested in solving problems of type ( [T4| under more general conditions, 
where the complete eigenvalue problem requires large computational costs. 

In this situation, we cannot directly apply the well-developed iterative methods 
of linear algebra and an appropriate software for solving ( p~4) >. This results from two 
reasons. On the one hand, we have the term pA°y on the left side. On the other hand, 
the equation is multi-term, i.e., it is represented as the sum of two individual operators. 

Obviously, elliptic problems with the operator A can be solved in an efficient way. 
Then the operator A can be selected as a preconditioner for iterative solving equation 
Let y k be an approximate solution at the k -th iteration. If we apply the conjugate 
gradient method El ED- then the new iteration is defined as follows. Denote r k = 
1 - Ay k , A - A + jjA" as the original residual and let z, k = A 1 r k be the residual for the 
preconditioned equation. With the initial po = zo and the given yo- for k - 0,1,..., we 
have 


Uk 


Zjfc+1 


fa, n) 

(Ap k ,pk) 

A~ l r k+ 1 , 


yk+ 1 =yk + a k p k , r k+ 1 =r k - a k Ap k , 


_ (Zk+\,r k +\) 
C Zk , n) 


Pk +1 - Zk+l + PkPk- 


(15) 


The convergence rate of the iterative method ( fl5] l is governed Il20l by the constants 
y\ and yz (more precisely, by the ratio k = 71 / 72 ) in the following bilateral operator 
inequality: 

71 A < A + pA a < 72 A, 71 > 0. (16) 


For A, in view of (12 1 and 0 < a < 1, we have 


A=A+pA a >A, A + pA a = {E + pA a ~ l )A < (1 + pd a ~ l )A. 


Therefore for 71 and 72 in ( fl 6 | , we obtain 

7i = l. 72 = 1 +p8*~ l . 


This establishes the dependence of the number of iterations in the conjugate gradient 
method ( fl 6 | on p, 6 and a. 

At each iteration, we must evaluate the quantity 

Ap k = A(p k + pA a ~ l p k ). 

The emphasis here is on calculating w = A" 1 p k . It is necessary to solve the problem 

A (i w = /, (17) 


where /3 = 1 — a, f — p k for 0 < p < 1. We apply the approach proposed in the paper 

ED. 

An approximate solution is sought as the solution of an auxiliary evolutionary prob¬ 
lem, where t is the pseudo-time evolution variable. Assume that 

v(t) = ( 96) a (t(A - 86E ) + 85Ey a v( 0), 
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with 0 < 9 < 1. Therefore 


v( 1) — (96) a A~ a v( 0) 

and then w = v(l). The function v(t) satisfies the evolutionary equation 

dv 


where 


By <[T2]», we get 


(tD + 961)—+ aDv = 0, 0 < t < 1, 

at 


D - A - 65E. 


D = D* > (1 - 6)6E > 0. 


We supplement equation (18 i with the initial condition 

v(o) = m~ a f. 


(18) 


(19) 

( 20 ) 


The solution of equation ( fT7] > can be defined as the solution of the Cauchy problem 
(p~8]>—(|20]> at the final time moment t = 1. In fl6l . the case of 9 - 1 was studied. 

For the solution of the problem ( p~8] >, ( |20| ), we can obtain various a priori estimates. 
Elementary estimates have the form 


IIvMIIg < llv(0)|| G , 


( 21 ) 


where, for instance, G — E,D. To obtain pT| for G — D, multiply scalarly equation 
( p~8] > by dv/dt. For G - E, equation ( p~8]> i s multiplied by av + tdv/dt. 

To solve numerically the problem (| 1 8|), ( |20| , we use a simple two-level scheme. Let 
r be a step of a uniform grid in time such that v" = v(t n ), f = m,n = 0, l,.... No, Not = 
1. Let us approximate equation ([T8| by the Crank-Nicolson scheme 


(t n+1/2 D + 66 E)- 


+ aD- 


= 0, n =0,1,..., No- 1, 


( 22 ) 


v° = m~ a f. 


(23) 


approximates the problem (18 1 , (20 1 with the second 


The difference scheme (22 1 , 
order by r. 

For 0 < 9 < 1, the difference scheme ( [22] ), ( [23] ) is unconditionally stable with 
respect to the initial data. The approximate solution satisfies the estimate 


l|v" +1 || G < ||v°|| c , n = 0,1, ...,No - 1, 


(24) 


with G — E,D. 


Multiplying scalarly equation (22 1 by v" +1 - v", we get 

l|v n+I |b < l|v"|b, n = 0,1, ...,Nq - 1. 


This inequality ensures the estimate (24 1 for G = D. 


Similarly, we consider the case with G - E. Rewrite equation (22 1 in the form 


96- 


v «+l _ v « 


+ D\a 


v n+ 1 + v" 


j.n+l/2 


v n+ 1 - v r ‘ 


= 0 . 
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Multiplying scalarly it by 


v" +1 + v n 


a- 


+ f 


n +1 /2 


v" + i - v" 


in view of (fi~9|, we arrive at 


< 0 . 


We have 

[|v" +1 || < ||v"||, n = 0, l,...,No - 1, 
that means the fulfilment of the estimate ([24} with G = E. 


4 . Numerical results 

To discuss our predictions, we start with calculations of the problem with the frac¬ 
tional power ( fl~7] >. The scheme (j22|, ( |23] > was applied. The problem, unless otherwise 
stated, was solved on the spatial grid N\ = Ni = 100, d = 1 with p = 0.5, 96 = 2n 2 . The 
evolution histories of the maximum value w max of the approximate solution (located at 
the center of the domain) are shown in Figure [T] for various computational grids in the 
pseudo-time evolution variable t (No - 5,10,20,100). It is easy to see that even on 
coarse grids in t, we observe a good accuracy of the solution. Figure [2] demonstrates 
similar pseudo-time histories calculated starting from other initial value of <[23l» in the 
scheme (22 1 . Namely, in this case, we use a rougher initial approximation 96 


7 T that 


corresponds to an inaccurate estimation for the lower bound ( p~2| ) of the operator A. 

The non-local convergence of the approximate solution with refining the grid in the 
pseudo-time evolution variable t is depicted in Figures [3] and J4lf or solution profiles. 
There are presented profiles of the solution of the problem ( |12| ) in the mid-section 
(X 2 = 0.5). As above, these profiles of w(xi,0.5) were calculated using various grids 
in t and starting with two different initial values. Obviously, the maximum error is 
observed in the vicinity of boundaries of the computational domain. 

The solution convergence for the fractional Laplace operator problem with refining 
the grid in space is shown in Figure [5]for the above mid-section profiles of w(xi,0.5). 
The calculations were performed on the finest grid in t (No = 100). 

For the problem ( fl7j ), the main interest is in the impact of the power p on the 
solution features. To eliminate the influence of grid parameters (grid steps in space 
and the pseudo-time evolution variable t), all predictions in this parametric study were 
performed on the finest grid with N\ — Ni = 100 and No = 100. Figure [6] presents 
mid-section profiles of the solution for various values of (3. For the convenience of a 
comparison, the solutions are normalized to the maximum value. It is easy to see that 
the decreasing of [3 leads to more gently sloping profiles. WhenyS —> 0, we have w —» / 
in the computational domain O. The dependence of the solution maximum w max on the 
power p and the geometry (the width of the computational domain d) is presented in 
Figure [7] 
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f the solution for/3 = 0.75 















Figure 11: The relative error s vs iteration number k for various values of fi (a = 0.5) 


The solution of the problem (17 1 normalized to the maximum value is shown in 
Figures [8|jT0| as isocontoures in the whole computational domain for different values 
of the power ft. We can observe the formation of a boundary layer when /? —> 0. 

Now we discuss the main object of our study, i.e., the problem ( [T4| . To solve 
it, we apply the iterative method of conjugate gradients with the operator A as 
a preconditioner. From the methodological point of view, the most interesting is the 
dependence of the iteration number on the parameters /j and a. The decreasing of the 
relative error s\ = ||r>||/||ro|| during iterations (with the initial approximation yo = 0 ) is 
given in Figure[TT]for various values of /j. The problem was solved with a = 0.5. The 
dependence of the convergence rate on a for /r = 100 is presented in Figure [l2| 

In modeling turbulent flows by means of the space-fractional model, we operate 
only with mean values of the longitudinal velocity component. A more detailed de¬ 
scription of turbulent flows is carried out on the basis of more complicated models of 
turbulence (see, e.g., EHE2). To validate our space-fractional model of turbulence, a 
comparison with experiments was done. A fully developed turbulent flow in a square 
duct was measured in |23l . We use experimental data from this study, which are placed 
on the Internet resource http://www.jsme.or.jp/ted/HTDB/fw.html Experimental pro¬ 
files of the normalized mean longitudinal velocity u mean are shown in Figure [13] for 
various cross-sections of x\ for a half of a cross-section. Here the origin of coordinates 
is located at the left bottom corner of the duct cross-section and so, at the center of the 
duct we have x\ = 0.5, X 2 = 0.5. We see more gently sloping profiles of the velocity 
in approaching to duct walls. Also we see increasing of the velocity towards the cor¬ 
ners of the duct, which is associated with secondary flows observed in experiments and 
which it is difficult to reproduce using simple models of turbulence. 

These experimental data we used to tune the parameters of our space-fractional 
model of turbulence (fl4|> in order to meet the above experimental data in the best way. 
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Figure 12: The relative error £ vs iteration number k for various values of a (ji = 100) 


For this purpose, a parametric study with respect to p and a was done. We estimated the 
deviation between the calculated and measured values of the longitudinal velocity. Let 
xl- 1,2,..., L be the points of measurement. The deviation measure is the quantity 


a) 


xi) - li(xi )) 1 


1/2 


V/=l 


where yip., a; X/) is the predicted velocity, whereas u(x/) is the measured velocity. Fig¬ 
ure 14 demonstrates the dependence of g on p for optimal values of a. These results 


show that the first term in the left-hand side of equation ( [T4| can be neglected. There¬ 
fore it is possible to use the one-term diffusion model, where instead of ||6}, we consider 
the equation 

dv 1 

— + v ■ Vv + - V/J + £(- A)“v = 0. 
at p 

For a flow in a duct, we can reduce equation ([8]) to the following equation 




(25) 


for the longitudinal velocity. 

Numerical results obtained using the one-term space-fractional model of turbulence 
( p5| at near optimal values of a are presented in Figures [ISj - fTTI The calculated data are 
compared with experimental profiles along three different lines of the duct xi = 0.5, 
xi = 0.7 and X\ = 0.9, respectively. A good agreement between approximate solutions 
and the measurements is observed in the critical region near duct walls. Relatively 
large discrepancies take place only in the central zone of the duct cross-section. This 
is partly due to the fact that the measurement points have a non-uniform distribution, 
i.e., near the boundaries the distance between the points is eight times lower than near 
the center. 
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Figure 13: Experimental profiles of u mean (23) along various cross-lines of x\ 



Figure 14: Dependence of ^ on n for optimal values of a 
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Figure 15: Comparison of experimental (solid line) and numerical profiles of u mean for x\ = 0.5 



Figure 16: Comparison of experimental (solid line) and numerical profiles of u mean for xi = 0.7 
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Figure 17: Comparison of experimental (solid line) and numerical profiles of u mean for x\ = 0.9 
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